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ABSTRACT 

Aims. The structure and potential of a complex gravitational lens is reconstructed using the perturbative method 
presented in Alard 2007, MNRAS, 382L, 58; Alard 2008, MNRAS, 388, 375. This lens is composed of 6 galaxies 
belonging to a small group. 

Methods. The lens inversion is reduced to the problem of reconstructing non-degenerate quantities: the 2 fields of 
the perturbative theory of strong gravitational lenses. Since in the perturbative theory the circular source solution is 
analytical, the general properties of the perturbative solution can be inferred directly from the data. As a consequence, 
the reconstruction of the perturbative fields is not affected by degeneracy, and finding the best solution is only a matter 
of numerical refinement. 

Results. The local shape of the potential and density of the lens are inferred from the perturbative solution, revealing 
the existence of an independent dark component that does not follow light. 

Conclusions. The most likely explanation is that the particular shape of the dark halo is due to the merging of cold dark 
matter halos. This is a new result illustrating the structure of dark halos at the scale of galaxies. 
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1. Introduction. 

In 1986 Lynds and Petrosian reported the discovery of gi- 
ant arcs in Abell 370 and 2242-02, with a length of 100 Kpc 
and luminosities comparable to giant eUiptical galaxies, 
but with bluer colors. Bohdan Paczyhski (1987) interpreted 
these features, he suggested that the arcs may be the result 
of the clusters acting as gravitational lenses. The Paczyhski 
interpretation was later confirmed by Grossman & Narayan 
(1988), and initiated a series of very fruitful developments 
on gravitational lenses. The lens inversion problem was 
tackled by Kochanek et al. (1989) using the ring cycle 
method. Other methods have developed to solve the lens 
inversion problem. Such as maximum entropy methods by 
Wallington et al. (1996), methods based on Fourier series 
expansion (Trotter et al. 2000), LensCLEAN (Wucknitz 
2004), semi linear inversion methods (Warren & Dye 2003), 
Bayesian methods (Brewer & Lewis 2006), and potential 
reconstruction method (Suyu et al. 2009). The lens inver- 
sion methods may be classified in 3 types: (i) model de- 
pendent reconstructions, (ii) potential reconstruction on a 
grid, (iii) expansion of the potential in Fourier series. All 
methods have to consider the degeneracy of modeling, for 
instance Wayth et al. (2005) found that 6 different mod- 
els were consistent with the observations for the lens ER 
0047-2808. The problem of degeneracy in lens inversion is 
widespread. The degeneracy issue is related to the fact that 
the constraints on the potential are local. The value of the 
potential is constrained only in areas where images of the 
source are formed. Additional constraints comes from the 
fact that no images are formed in dark areas (Diego et 
al. 2005). Since there are no constraints on the potential 
in other areas, without additional assumptions, the prob- 



lem is fully degenerate, and an infinite number of mod- 
els are possible. One possibility is to explore the family of 
models that is consistent with the observations (Saha & 
Read, 2009) and to try to find the common features be- 
tween the models (non-degenerate quantities). However, in 
such an approach it is difficult to explore the full range 
of possible models. Another possibility to reduce the de- 
generacy is to use a local representation of the potential, 
like the methods (ii) and (iii). However, the grid methods 
or Fourier methods require many free parameters, and it 
is not clear whether the degeneracy problem has been im- 
proved. To solve the degeneracy problem it is essential to 
reduce the number of parameters, and to describe the po- 
tential in terms of fundamental, non-degenerate quantities. 
Since giant arcs resemble a part of a circle, it is natural 
to consider that arcs corresponds to a small break of the 
circular symmetry. In such an hypothesis, the fundamental 
properties of the arcs should be related to the quantities 
breaking the circular symmetry. These ideas were first ex- 
plored by Blandford & Kovner (1988) who demonstrated 
that for a small impact parameter (with respect to the size 
of the critical circle) and small deviations from circular- 
ity, the equation for caustics depends only on the orthora- 
dial component of the defiection field at the critical circle. 
Blandford & Kovner (1988) also presented a geometrical 
method to reconstruct the images of the lens. The same con- 
cepts were investigated later by Alard (2007) who derived 
an analytical equation to describe the formation of images 
by the lens. The main advantage of the analytical descrip- 
tion is that it directly relates the morphological properties 
of the arcs to the potential expansion. This relation between 
the potential and the observational features suggests that 
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the description of the potential is non-degenerate, and also 
simplifies the lens inversion problem. The reduction of the 
degeneracy problem in this approach is well illustrated by 
the fact that an infinite number of models corresponds to 
a given perturbative representation. Lens inversion in this 
perturbative framework was already tested in Alard (2008) 
using numerical simulations. It was found that the pertur- 
bative approach allows an accurate re-construction of the 
potential for lenses with substructures. The concepts devel- 
oped in this approach are very different from some recent 
perturbative numerical schemes (see for instance Vegetti 
& Koopmans 2009). In this type of work there is no at- 
tempt at reformulating the lens-inversion problem in terms 
of fundamental quantities; the main goal is to derive an effi- 
cient numerical scheme to solve the ordinary lens equation. 
Peirani et al. (2008) showed that the perturbative approx- 
imation is accurate for realistic lens models. The numeri- 
cal accuracy is of about 1% in units of the critical radius. 
The lens model of Peirani et al. (2008) were extracted from 
cold dark matter simulations (Horizon project), and some of 
them have a lot of structure (merging) . This paper is orga- 
nized as follows: the first part will present the perturbative 
method, and its basic properties. In an effort to make this 
work self-contained, a synthesis of the other useful results 
on the pertubative method is presented in the appendices. 
This section also has the advantage of presenting all the 
equations essential to lens inversion. The second part of 
the paper will present the application of the perturbative 
method to the inversion of an interesting case of gravita- 
tional lens. This paper also will be an occasion to extend 
the technique of perturbative inversion and to relate the 
perturbative fields to the geometry of the lens. 

2. The perturbative approach in gravitational 
lensing. 

In the perturbative approach (Blandford & Kovner 1988 
and Alard 2007), gravitational arcs are interpreted as a 
perturbation of the perfect ring situation. A perfect ring 
is formed when a point source is situated at the center of 
a circularly symmetric lens. In the perfect ring situation, 
the central point has an infinite number of images located 
on the critical circle. There are 2 different kinds of per- 
turbations that breaks the symmetry of the system: the 
non-circularity of the lens, and the source impact parame- 
ter. When such non-circular perturbations are introduced, 
the symmetry of the system is lost and the ring is replaced 
with a series of individual images. However, as illustrated 
in Alard (2007) whatever the angular position 6 of these 
perturbed images, there is always an unperturbed image of 
the central point on the critical circle at the same angular 
position. As a consequence there is only a radial displace- 
ment dr between the perturbed and un-perturbed images. 
The un-perturbed situation is represented by a circular lens 
with potential 0o('')j and a point source with null impact 
parameter. In the perturbed situation, the source has an 
impact parameter rs , and the lens is perturbed by the non- 
circular potential ^(r, 9). Both perturbations rs and ip are 
assumed to be of the same order e: 



(1) 



Note that for convenience the unit of the coordinate is the 
critical radius. Thus by definition in this coordinate system 



the critical radius is situated at ?' 1. As a consequence, in 
the continuation of this work all distances and their associ- 
ated quantities (errors,..) will be expressed in units of the 
critical radius. In response to this perturbation, the radial 
position of the image will be shifted by an amount dr with 
respect to the unperturbed point on the critical circle. The 
new radial position of the image is: 



r = 1 + e dr 



(2) 



The response dr to the perturbation can be evaluated by 
solving the lens equation in the perturbative regime. The 
lens equation reads: 



Ts = 



dr 



r'de 



(3) 



The last step to solve the lens equation is to expand the 
potential is series of e: 



- ^ [C„ + e Um [r - 1) 



n=0 



with: 




(4) 



(5) 



By inserting Eqs ([T]), ([J), and (2]) in the lens equation (Eq. 
[11) , it is straightforward to obtain an equation relating the 
response dr to the perturbation rg and "0: 



(K2 dr - fi) Ur - — U0 



(6) 



and K2 = 1 — 



(r=l) 



This equation corresponds to Eq. (8) in Alard (2007). 

Considering that the source has a mean impact param- 
eter To , the position in the source plane may be re- written: 
rs = rs +ro. Assuming Cartesian coordinates {xo,yo) for 
the vector ro, Eq. [6] reads: 



(^2 dr ~ /i^ Ur 



dfo 

de 



(7) 



For a circular source with radius _Roi the perturbative re- 
sponse dr takes the simple following form (Alard 2007, Eq. 
12): 



dr — — 

K2 



./i± 



7?2 - 



dfo 

de 



(8) 



Eq. ( [5]) depends on K2. However this variable can be 
eliminated from Eq. ( [6]) by re-normalizing the fields: /„ = 
— , and the source plane coordinates, r^ = — (mass sheet 
degeneracy) . These re- normalized variables will be adopted 
in the rest of this work. The re-normalization is equivalent 
to solving Eqs ( ^ and ( [5]) for K2 = I. The variable K2 will 
re-appear when the re-normalized quantities are converted 
to the original quantities. 
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3. Pre-processing of the HST data. 

The gravitational lens SL2SJ02140-0532 (Cabanac et al. 
2007, SL2S public domain) was observed by HST in 3 
spectral domains, F475W, F606W, F814W, with an expo- 
sure time of 400 seconds. The main difficulty in the pre- 
processing of the images is to remove the large number 
of cosmic present in the images. Cosmic are identified us- 
ing a wavelet approach based on an estimation of the local 
scale in the image. Once a cosmic has been identified at a 
given position, the local area around this position is inves- 
tigated to remove any other pixels connected to the cosmic. 
This procedure is based on routines from the ISIS package 
(Alard 2000). Note that since the number of cosmic is large 
in the images, the local density of cosmic sometime become 
too large in some areas. When too many nearby cosmic 
are detected within a small area, the area is flagged, and 
will be considered cautiously, or even rejected. Once the 
3 images are cosmic cleaned, a re-centering to a common 
grid is performed using the position of the bright objects. 
The last step in the processing of the images is to remove 
the background. There is some diffuse light coming from 
the galaxy halos at the level of the arcs, and the proper 
removal of this diffuse light is especially important. The 
diffuse background is removed in several steps. The first 
step is to build a general background model. To build this 
model, the image is divided in square areas of 32 pixel, the 
local background is estimated within each area using a kct 
cleaning method. The background value in each local area 
is then interpolated to compute a general background value 
at each position in the image. Once this background model 
has been subtracted from the image, the residuals around 
the arcs are investigated in order to control any systematic 
errors. A tiny systematic background residual is detected; 
this residual depends on the distance to the 3 main galaxies 
center of gravity. This background residual is related to an 
imperfect subtraction of a diffuse halo centered around the 
3 main galaxies. An elliptical component centered at the 
center of gravity of the 3 galaxies is fitted to this residual. 
The best elliptical parameters are estimated by minimiz- 
ing the deviations to elliptical symmetry. The mean radial 
profile of the residual is estimated, and finally the relevant 
elliptical component is subtracted. Note that this last back- 
ground correction is small, it is only a few percent of the 
initial background subtraction. Thus, the 3 cosmic cleaned, 
re-centered and background subtracted images are stacked, 
to produce the final reference image of the arc system. A 
color illustration of the arc system is provided in Fig. [T] 

4. Lens inversion. 

This work presents the first successful model of the images 
formed by the lens SL2SJG214G-G532. This system was in- 
vestigated by Limousin et al. (2009), who tried to fit this 
system using an elliptical lens model. They found that it 
was not possible to reproduce this system using such simple 
elliptical models. Limousin et al. (2009) were also confused 
by the color of the last image (right side in Fig. [1]) which 
is slightly different to the color of the other images. They 
concluded that the identification of the images formed by 
this lens was uncertain, and as a consequence renounced 
modelling of this system. We will see in the continuation of 
this work that this lens requires models much more complex 
than elliptical, and that such models can be built naturally 



Fig. 1. Color image for SL2SJ02140-0532 reconstructed 
from 3 noise filtered HST images in 3 bands. 



using the perturbative method presented in this paper. The 
development of this perturbative model also will naturally 
solve the problem of the different color of the last image. 
We will see that the first part of this image has the same 
color as the other images, while the remaining part of the 
image has a different color. This is a direct prediction of 
the perturbative model, which offers an elegant solution to 
the problems encountered by Limousin et al. (2009). We 
now turn to the analysis of this lens using the perturbative 
method. The lens SL2SJ021408-053532 forms a system of 4 
images of the source. The first 2 images are situated inside 
the large arc in the upper left (they appear as 2 identical 
mirror images). The 2 other image, are at the bottom and 
the right side of the frame (see Fig [T|) . In the perturbative 
approach the problem of reconstructing this lens is equiv- 
alent to the reconstruction of the perturbative fields, and 
the next sections will concentrate on this issue. Obviously, 
the re-construction problem is reduced to the reconstruc- 
tion of the first order perturbative field only if the first 
order perturbative theory applies. The requirement for this 
approximation to work is that the gradient of the potential 
can be linearized in the vicinity of the critical circle (Alard 
2007). For single dark matter halos numerical experiments 
shows that this local constraint on the gradient is usually 
met with good accuracy (Alard 2007, Peirani et al. 2008). 
It is only in the case where a local minima of the poten- 
tial is present close to the critical circle that this approx- 
imation may not work. But such situations are very rare 
because they require that for instance a sub-structure falls 
very close to the critical line. However Alard (2008) showed 
that the general treatment of substructure in this approxi- 
mation is correct even when the substructure is quite close 
to the critical circle (up to one tenth of the critical circle). 
Additionally, when an unlikely situation happens, such as 
a substructure lying on the critical line, the consequences 
are easily observable: a few additional images are formed 



4 



C. Alard: Perturbative reconstruction of a gravitational lens: when mass does not follow light 



near the perturbator. Thus such situation is easy to recog- 
nize, and is not apparent in the present data. In any case 
the best test of the first order approximation is to perform 
a successful reconstruction of the data, which means that 
in practise second order terms can be neglected. This will 
be the approach favored in this paper. Before the field re- 
construction problem can be tackled, a necessary step is to 
evaluate the size of the critical circle, which will become by 
convention the coordinate system unit. 

4.1. Estimation of the critical circle. 

The critical circle is estimated by fitting a circle to the mean 
position of the images. The center and radius of the circle 
are adjusted. The non-linear adjustment procedure starts 
from a circle centered on the small galaxy at the center 
of the image. The initial estimate for the circle radius is 
the mean distance of the images to the center. A non linear 
optimization shows that the best center is close to the initial 
guess and that the optimal radius is close to 152 pixels; with 
a pixel size of 0.049" the critical circle radius Rc is 7.44". 
This value should be compared to the radius of 7.31" found 
by Limousin et al. for the same lens. Note that the final 
result will not depend upon the particular choice of a given 
critical circle. Taking another circle close to this one would 
change a little the estimation of the perturbative fields, but 
the total background plus perturbation would remain the 
same. 

4.2. General properties of the solution. 

The general properties of the solution can be derived from 
the properties of the circular source solution. To estimate 
the circular source solution, the outer contour of the source 
will be approximated with a circle. This approximation of 
the source may not be very accurate; the typical error on 
the field estimation will be of the order of the source outer 
contour deviations from circularity. However, this approxi- 
mation is sufficient to obtain a first estimate of the solution 
and to derive its general properties. For circular sources 
the perturbative fields are directly related to the data (see 
Appendices El and Eq. [8]). This direct relation between 
the fields and the data is possible only in areas of the lens 
plane where images of the source are present. In dark areas, 
the field will have to be reconstructed by interpolating the 
values in bright areas. The field /i will be approximated 
with the mean radial position of the images (with respect 
to the critical circle). The reconstruction of /i is simple and 
not ambiguous, but the reconstruction of ^ requires more 
work. On the other hand ^ is better constrained since 

do 

no images should be formed in dark areas (|^| > ^o)- 

Furthermore, the field is of particular interest since it 
is directly related to the potential iso-contour (Eq. ( IC.lj) ). 

The derivation of the general properties of ^ require some 
specific guess of the local behavior of this field in the region 

of image formation. Images form in minima of 1^1 and the 
example presented in Appendix |X] indicates that the be- 
havior of ^ near the minima is of 2 types. For small images 

has a linear behavior, while for larger images (caustics) 
the field behaves like a higher order polynomial. In the case 



of this particular lens system, we have 4 images of compa- 
rable size. The typical image size is quite small and it is 
reasonable to assume that the field behavior in the image is 
linear. Of course other models could be investigated, with 
for instance order 2 behavior of the field, but obviously 
these models are much more complex and would require 
higher order Fourier expansions (see Fig. [2]). Additionally, 
adopting a higher order description (lower curve in Fig. [5]) 
may introduce false features like local peaks. The locally 
linear model (upper curve in Fig. [2]) is the simplest de- 
scription consistent with the general properties of ^; it 
has the smallest possible number of minima and maxima. 
Thus adopting the locally linear model guarantees that no 
un-necessary features are introduced in the model. On the 
other hand, in reality this description may be too simple. 
The validity of this model can only be tested by a direct 
application to the data, and an evaluation of the quality of 
the fit. A higher order description should be adopted only 
if the locally linear model does not fit the data properly. 
As a consequence, the locally linear model will be adopted 
in the continuation of this work. Once the general shape of 
the solution is known, everything is now just a matter of 
improving the numerical accuracy. The process of building 
the solution will proceed in several steps. First the local so- 
lution will be refined numerically, then some interpolation 
of the solution in dark areas will be computed. A Fourier 
expansion will be fitted to this piecewise polynomial solu- 
tion. And finally the adjustment of the Fourier expansion to 
the data will be optimized using a non linear minimization 

procedure. Note that since only the square of ^ enters 

Eq. ( [8]), the sign of ^ is not known, as a consequence 
two solutions with different sign will have to be explored. 
Hopefully, when the source is not circular this degeneracy 
of the sign can be broken. 




Fig. 2. Two example of possible shapes for the field -j^. 
The black dots represents the mean image position. The 
red curve represents the local behavior of the field near 
the image, while the black curve is a smooth extrapolation. 
The upper curve present a solution where the field has local 
linear behavior near the images. In the lower curve a higher 
order behavior is presented. In dark areas the field must be 
larger that the source radius to avoid forming images (see 

Eq. [8]). Note also that the field ^ must also have zero 
sum. 
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4.3. Refinement of the local solution. 

The study of the general properties of ^ indicates that the 
simplest model of this field is locally linear in the neighbor- 
hood of the images. This local linear model will be refined 
in this section. The same local linear approximation will be 
used for the field /i . As illustrated in the former section, the 
estimation will start from the circular source guess. For a 

circular source and a linear field model, the field ^ is zero 
at the center of the image. This property will be assumed 
in the refinement of the local solution. The slope of the lin- 
ear model for ^ is estimated by solving Eq. ( [8]). Note 
that solving Eq. ( [H]) requires an estimation of the source 
radius Rq. Since the image maximum half radial width is 
i?0) an estimation of Rq can be obtained by averaging the 
maximum radial width of the 4 images. By adopting this 

linear local model of taking the mean radial position 
as an estimation of /i linear parameters, and using the per- 
turbative lens equation (Eq. [H]), it is possible to fold the 
4 images to the source plane. Once folded to the source 
plane, the 4 images must be similar and their comparison 
is a direct test of the quality of the lens re-construction. 
With the initial guess for the field a, correlation between 
the 4 folded images is effectively observed: the mean cross- 
correlation between the images is 0.6. This cross correlation 
can be improved quite significantly by non-linear optimiza- 
tion. The cross-correlation between images in the source 
plane is maximized using the Nelder & Mead (1965) sim- 
plex method. The parameter space of the maximization is 

the 4 local slopes of the ^ field, and also the 4 local slopes 

of the field /i . The simplex maximization brings the mean 
cross-correlation between images to 0.8 which is a good in- 
dication that the linear approximation is already a good 
model. This maximization is also the occasion to build a 
first realistic model of the source, and to go beyond the 

circular source approximation. The result for the field 
is presented in Fig. ( [3]) . The local linear approximation is 
the dark blue line segments in Fig. ([3]). 

4.4. Interpolating the local linear solution in dark areas. 

Since there are no constraints on the field /i in dark areas, 
the fit of this field in the full angular range is straightfor- 
ward. Due to the constraint that no images are formed in 
dark areas, the situation is more complicated for the field 

The local linear solution obtained in the former section 
needs to to be completed with second order polynomials to 
fill the gaps between the images. The condition that no im- 
ages are formed in dark areas (Diego et al. 2005) is equiv- 
alent to the condition that the functional |^| is larger 
in dark areas than in areas of image formation. Another 

constraint is that the sum of the field H§- must be zero. 

do 

And finally the field must be continuous and as smooth as 
possible. The first constraint is a constraint on the sign of 
the curvature of the local polynomial, or equivalently on 
the sign of the second order coefficient of the polynomial. 
The second constraint reduces the number of unknowns in 
the polynomial coefficients. And the last constraint settles 
the value of the remaining coefficients. Note that the con- 
straints on zero and first order polynomial coefficients can 
be solved using linear methods, which reduce the number 
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Fig. 3. An estimate of the ^ field. The blue hne segments 
represent the local linear guess at the image location, while 
the light blue curves are the second order interpolations 
between images. The red line is an adjustment of a 6th 
order Fourier series to the blue curve. 



of unknown and leaves us with variables depending only on 
the second order coefficients. The remaining equations de- 
pending on the second order coefficients have to be solved 
with a constraint on the sign of the coefficients. A simple 
grid method was adopted to solve the relevant equations. 
The final solution is presented in Fig. [3] (light blue curve). 

4.5. Final fitting of the fields. 

The final fitting of the fields will be conducted by recon- 
structing the images of the source and comparing these im- 
ages with the HST data by chi-square estimation. This pro- 
cedure will take into account the convolution of the images 
with the PSF, which was ignored in the former reconstruc- 
tions. The HST PSF is estimated using the Tiny Tim soft- 
ware (Krist 1995). Note that the chi-square value between 
model and data can be interpreted as a comparison between 
the images. In the model, all images come from the same 
source, thus the chi-square directly measures the image sim- 
ilarity in the lens plane. In practice, the fitting procedure 
will work in the following way: start from a numerical guess 
of the field, reconstruct the source, then reconstruct the im- 
ages of the source, and finally estimate the chi-square with 
respect to the HST data. An essential component in this 
procedure is the source reconstruction method, which is il- 
lustrated in the next section. 

4.5.1. Source reconstruction 

At this point the simplest method to estimate the source 
would be to use the Warren & Dye (2003) method. 
Unfortunately the source is very complex, and this method 
would require using too many basis functions, which in 
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practice is un-tractable. The following method has been 
preferred: images are co-added in the source plane using an 
algorithm that preserves the geometrical properties of the 
images elements. The basic idea of this geometrical method 
is that pixels in the lens plane are divided in 2 triangles. 
These triangles are transported to the source plane and are 
assumed to remain triangles. The triangles in the source 
plane are projected on a grid and the flux is co-added. The 
optimal size of the grid is estimated by numerical experi- 
ments. The noise level in the source image is estimated by 
looking at a specific area, and only the significant pixels 
in the source grid are kept. The images of the source in 
the lens plane are computed using the same geometrical 
method, but this time in the opposite direction. The final 
step is to compare the images of the source with the data. 
The drawback of this method is that the pixels are affected 
by the convolution with the PSF, and thus the source re- 
construction is not perfect. However, it is easy to correct 
this effect by using the following procedure: first ignore the 
effect of the PSF on source reconstruction, find a solution, 
and then correct the effect using a deconvolution procedure 
very similar to the Van-Cittert method (Van Cittert 1931). 
This de-convolution of the source starts from the solution 
obtained ignoring the PSF convolution. The initial guess is 
convolved with the PSF and a difference image is formed 
by subtracting the data. This difference image is taken as 
a new input to the source re-construction procedure in or- 
der to produce a first correction of the initial guess of the 
source. This procedure is further iterated until only noise 
is present in the difference image. In practice this method 
requires to control the noise in the difference image. In this 
work a 3-a level was adopted; tests were performed with 
other thresholds, 2 and 4 a, with no significant change in the 
results. The lower a cuts require more iterations and pro- 
vide somewhat more noisy source re-construction. However, 
since the areas where effective corrections are needed are 
quite small (these areas correspond to the bright sharp fea- 
tures in the images), the source reconstruction is not af- 
fected much. It is also important to note that these PSF 
corrections have very little effect on the re-construction of 
the perturbative fields. By comparing the fit obtained by 
ignoring the PSF and correcting for the PSF, it was found 
that the difference is only a small fraction of the noise in 
the field re-construction. Thus for this particular lens this 
procedure of PSF de-convolution is important only for the 
re-construction of the source. This is due to the fact that for 
this particular lens the thickness of the arc system is quite 
large with respect to the PSF size, and as a consequence 
the PSF convolution has only a limited effect. 

4.5.2. Fitting of a Fourier model. 

The practical implementation of the fitting procedure pre- 
sented above now will be described. The first step is to 
choose a starting estimate. It could be the piecewise poly- 
nomial model presented in Sec. 4.6, but since the Fourier 
expansion of the fields is directly related to the multipole 
expansion of the potential (see Eqs [RT] |R2]and |R3)) . it 
is more convenient to use Fourier series. The initial guess is 
obtained by fitting a 6th order Fourier series to the piece- 
wise model of and for /i by fitting a Fourier series to 
the mean radial position of the images. Starting from this 
initial guess, the source and image reconstruction proce- 



dures are iterated using the simplex method, in order to 
minimize the chi-square between the model and data. This 
procedure converges to the minimum of the chi-square in 
the parameter space. This minimum corresponds to the op- 
timal Fourier coefficients. Note that as pointed out in Sec. 
4.2, the sign of the initial guess is unknown and 2 solutions 
with different signs have been explored. Fortunately in this 
case the source is not circular and the 2 solutions have 
different chi-square, which solves the problem of the degen- 
eracy of the solution sign. The final solution, the source 
reconstruction and the image reconstruction are presented 
in Fig [5l and the field solution is presented in Fig. (|4]). The 
fact that the source is complex and has a number of small 
features (see Fig. [7]) can be used to estimate the accuracy 
of the modeling. There are 3 bright features in the source 
that are visible in all images (red, yellow, blue in Figs [6] 
and [7]). The position of these features in the HST data and 
model has been compared for each image. To obtain the 
best possible accuracy in the comparison of the positions, 
the following procedure was used: a small image of each fea- 
ture is extracted in the model frame, and this small image 
is cross-correlated with the HST data at various positions 
around the feature. The cross-correlation map obtained us- 
ing this procedure is fitted with a polynomial to locate the 
map maxima, to derive a direct estimate of the shift be- 
tween the features. The residuals are somewhat smaller for 
the first 2 images. In particular a small discrepancy is visi- 
ble in the bright part of the third image (lower left image); 
the distance between the two images of the bright features 
in the source are closer in the model. The discrepancy is of 
about 2.4 pixels, which is about 1.6 % of the critical radius. 
On average the mismatch between model and observations 
is 0.9% of Rc- This estimation is very close to the results 
obtained by Peirani et al. (2008) who found that for re- 
alistic lens models the mean errors on the reconstruction 
of images using the perturbative approximation was about 
1% of Rc- Note that two other features are visible in the 
source (see Fig. [7]); due to variable levels of amplification, 
noise, and blurring by other nearby features, the images 
of these features are not always clearly visible in the HST 
data. However, when they are visible in the images, their 
position and the correspondence between model and data 
are indicated in Fig [S] Note that the caustic system asso- 
ciated with the solution is complex and includes a number 
of sharp features (see Fig. [5]). The complexity of this caus- 
tic system is not easily visible in the images; this is due to 
the fact that the source is large with respect to the small 
features in the caustic system. 



4.6. Noise. 

The model and the data are compared in a small region 
around the arcs with total number of pixels N. The Poisson 
weighted difference Ri between the model at pixel i, Mi 
and the HST data Di is very close to a Gaussian (see Fig. 
[8]). Considering that the model has Np parameters, 25 pa- 
rameters of the Fourier expansion of the fields, and the 
approximately 450 pixels of the source, the chi-square is, 
X2/dof = N-Np S» - l-l'''- Changing a little the size of 
the area, either reducing it by moving closer to the center of 
the arcs, or enlarging the area does not significantly change 
the chi-square value. 
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Fig. 5. Results from image re-construction: HST data (upper right), model (upper left), data and model subtraction 
(lower left), source (lower-right). Note the reproduction of the source details in the images of the model. An image of the 
source has been super- imposed on the model (upper left), to illustrate the scaling. 



4.6.1. Errors on the Fourier coefFicients due to noise in the 
data. 

The errors on the estimation of the Fourier coefficients are 
evaluated by local linearization methods near the minima 
(see, for instance Press et al. 2007). The intensity Mi is a 
function of the Np Fourier parameters p„, Mi = Mi{pn). 
Assuming that p„ is very close to the solution and adding 
a small shift to the parameters dp„, the intensity variation 
can be linearized. 



dMi 



> - — dp„ 



(9) 



The parameters dpn can be estimated by a linear least- 
square fit from Eq. ( (9]) , and thus the relevant errors are 
also estimated by taking the diagonal elements of ^ = 
, the normal least-square matrix invert B, Bnm — 



dMj 



This error estimation indicates that 

the errors on the Fourier coefficients due to the Poisson 
noise are ~ 0.5 10""^. 



4.6.2. Errors on the reconstruction of the perturbative fields. 

In the former section it was shown that the accuracy in the 
reconstruction of the images was close to 0.9 % of Rc- This 
error is directly related to the errors on the reconstruction 
of the perturbative fields. To estimate the amplitude of the 
errors on the fields it is particularly useful to consider cir- 
cular sources. For a circular source there is a direct relation 
between the errors on the fields and the errors on image 
reconstruction. Eq. ( [8]) shows that for circular sources, /i 
is the mean radial position of the image. As a consequence, 
the errors on /i are of the same order as the errors on the 
image position, which were estimated to be of the order of 

a percent. The field ^ is related by Eq. ([5]) to the radial 
width W of the image: 



Rl 



dfo 
d9 
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Fig. 6. Correspondence between the features in the model (right side) and the HST data (left side). The color code 
of the relevant source feature is presented in Fig. [T] The green ellipses indicate areas where cosmic cleaning was not 
satisfactory. There are a very large number of cosmic in the images, and sometimes several of them falls nearly in the 
same area, making cosmic cleaning very difficult. 



0.3 




1 2 3 4 5 6 

e 



Fig. 4. Best solution for the fields: red curve, and black 
curve"/!. The blue dashes indicates areas where images are 
present. 

As a consequence, the error on W and the error on ^ are 
related: 




Assuming the statistics of the variable W and are the 
same and differ only by a scale factor, the above equation 
shows that this scale factor must be unity. As a conse- 
quence, the errors on ^ are of the same order as the er- 
rors on the image radial width W , which are similar to the 
errors on image reconstruction: ~ 1%. The corresponding 
errors on the Fourier coefficients of the field expansion of 
order m can be estimated by analyzing the variance of the 
expansion: 

m 

Fi = ^ flj cos{j9) + bj s'm{j9) 
i=i 

Note that here Fq = ^ and Fi = /i . Taking the square of 
Fi and averaging over 9 




Provided it is assumed that the variance of (a|) = (b'j) = 
a^, the variance of Fi reads 

(F^) = ma^ ~ 0.01^ (10) 

Note that if the Fourier expansion has no constant term, the 
mean value of Fi is zero (/ Fid9 = 0). In practice there is an 
additional constant term for the field /i, however consid- 
ering that the number of Fourier parameters is quite large 
(12 parameters), the variance due to this additional param- 
eter will be neglected. As a consequence for both fields the 
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errors on the Fourier coefficients is ct ~ 0.004. This error is 
about 10 times larger than the error due to the Poissonian 
noise in the image, which is clearly negligible. As a conse- 
quence, in the continuation of this work, only the noise due 
to the perturbative approximation will be considered. 

4.6.3. Effects of the errors on the potential and density 
estimation. 

Some interesting properties of the lens are directly related 
to the perturbative fields as illustrated in Eqs IC.ll and 
ID. 21 The potential iso-contour depends on /o, and it will be 
assumed that the variance of /o is identical to the variance 
of /o. The variance on /o can be estimated using the method 
presented in Sec. 14.6.21 The expansion for /o read: 

-^^^^i cos(j6i) + bj ain{je) 

Thus 

/o = ^^cos(j0) + ^sinO-0) 
Using Eq ( [T0| the associated variance of /o is 

As a consequence, the average error on the estimation of /o 
is 5(piso — 0.005. In the same spirit the errors on the density 
iso-contours Spiso (Eq. ID.2P are 

, 2 0.01^ A ,,2 -2 

In practice the expansion in the above formula is conducted 
to order m, with 1 < to < m in order to limit the amplitude 




Normalized deviation 



Fig. 8. Histogram of the difference image between the 
model and data. The black curve is the histogram of the 
pixels in the difference image. For each pixel the difference 
has been weighted by the local Poisson noise expectation. 
The red curve is the theoretical Gaussian expectation for 
the normalized residuals. The closeness between the black 
and red curves indicates that the result of the fit is very 
close to optimal. The blue curve is the histogram of normal- 
ized residuals of the data without subtracting the model. 
This curve is plotted here for reference; it illustrates the 
effect of subtracting the model from the data. 



of the noise. Some numerical experiments shows that for 
TO > 4 there is no significant changes to the iso-contour 
shape, while on the other hand the noise grows quickly with 
TO. Thus taking m = 4, the resulting error is bpiso — 0.077. 




Fig. 7. The source features visible in the images (see Fig. 
[7]). The 3 brightest features are outlined by red, yellow, 
and blue ellipses. These features are used to estimate the 
accuracy of the reconstruction. The part of the source at 
the lower left of the image is effectively visible only in one 
of the images (see Sec. 5 for more explanations). 



5. Prediction of image color. 

Figure [§] indicates that the last image (lower right position 
in Fig. [T]) does not have the typical "hat-shaped" color 
diagram observed for the 3 other images. But this feature is 
actually a prediction of the perturbative model. The model 
predicts that the first half of the last image comes from the 
same source area as the other images. The other half of the 
last image corresponds to an area of the source that has no 
counterpart in the large arcs, and only a tiny counterpart in 
the lower left image. This situation is due to the crossing of 
the caustic line by a faint part of the source which is visible 
in Fig. [9] in the lower left of source-caustic diagram. The 
part of the source situated inside the caustics corresponds 
to the first half of the last image, while the other half of this 
image corresponds to the faint part of the source outside 
the caustics. As a consequence, the first half of the last 
image should have the same color diagram as the other 
images, and this is exactly what is observed (see Fig. [9]). 
The remaining part should have the same color as the faint 
parts of the source, and this is also what Fig. [§] illustrates. 
The color variations in the source are entirely related to the 
different colors of the bright and faint parts as illustrated 
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Fig. 9. The 3 color diagrams obtained by taking the sequence of raw flux in the 3 HST bands, as well as the source- 
caustic system (lower right). The upper left presents the color digram of the last image (red), and of the 3 other images. 
The upper right presents the color diagram for the bright parts of the first 3 images (red curves), while the blue curves 
represent the faint parts of the source. The lower left plot, compares the color diagram of the first half of the last image 
(black curve), to the mean diagram of the 3 other images (blue curve), and also compares the color diagram of the second 
half of the last image (green curve), to the mean color diagram of faint source parts in the 3 other images (red curve). 
The last figure presents the source and the system of caustics in the source plane. 



in Fig [HI The accurate prediction of the observed image 
colors is a confirmation of the model. 

6. Shape of the dark component 

The shape of the potential iso-contour is given by Eq. IC.ll 
However, in practice, only /o is known, and Eq. ( [7]) shows 
that the first order coefficients of fo depend also on the im- 
pact parameters. But these first order terms are related to 
the centering of the potential, not to its shape, for instance 
the ellipticity is related to order 2 Fourier terms. It is how- 
ever possible to obtain some partial information about the 
centering of the potential by using Eq. ( IB.3|) the first order 
moments within the unit circle do not depend on the im- 
pact parameters (they cancel out in the calculations). The 
numerical value of the inner moments is quite small, of the 
order of O.lRc- The outer moments have to be of the same 
order, and thus the centering terms are not large. Note that 
equations have been normalized by the factor K2 which de- 
pends on specific properties of the background. In practice 
K2 acts as an unknown scaling factor on fo and thus on 



the potential iso-contours. However, an important point is 
that this scaling degeneracy does not affect the shape of 
the potential, and for realistic halos K2 is close to unity 
(Peirani et al. 2008). The potential iso-contours dr ~ —fo 
are presented in Fig. f [TU|) . 

6.1. Mass to light relation 

The un-perturbed potential iso-contour is a circle with ra- 
dius r = I when a local perturbation is introduced, some 
distortions to the circular iso-contours will be visible in the 
neighborhood of the perturbator. As a consequence, pro- 
vided that the potential is generated by the visible matter, 
the distortions should correlate with the position of the 
galaxies. Since the center of the potential is unknown, the 
center of the un-perturbed unit circle is also unknown, and 
will have to be adjusted. The adjustment of this local circle 
is performed in the neighborhood of points on the potential 
iso-contour. The angular range of the circle adjustment is 
■|. The residual of this local unit circle adjustment is pre- 
sented in Fig. [ini The local deviations from circularity are 
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not consistent with the positions of the galaxies. The oppo- 
site is visible; the maximum deviations from circularity are 
observed in the lower part of the figure where no nearby 
galaxies are present. The average deviation from circular- 
ity of the potential in this area is ~ 0.05 Rc- According 
to the mean error on the potential estimated in Sec. 4.6.3, 
this is a lOcr deviation. This result demonstrates the need 
for a dark component that does not follow light. This is an 
interesting point for theories that try to avoid dark matter 
by modifying gravity. A natural solution is to consider that 
the dark component was formed by merging of cold dark 
matter halos. 

6.2. Inner and outer contribution to the potential. 

Note also that Eq. IB. 31 relates the Fourier coefficients of 
the fields ^ and fi to the inner (a„,5„) and outer terms 
{cn,dn) of the multipole expansion. As a consequence, the 
expansion of /o can be separated into inner contributions 
(which is equivalent to taking c„ = rf„ = in the original 
expansion), and outer terms (a„ = 6„ = 0). A numerical 
estimation shows that about 75 % of the total power spec- 
trum of the field /o comes from inner terms, which demon- 
strates clearly that the observed effects are dominated by 
the inner contributions. Since the deviation from circularity 
of the potential in areas with no light are much above the 
noise (lOcr), the result that mass does not follow cannot be 
affected by the 25 % contribution of the outer terms. 



some degeneracy (Clowe et al. 2006). The galaxies in this 
small group have very much the same colors, and as a con- 
sequence are probably very nearly at the same distance. 
Interpretating these results is not ambiguous and we are 
probably observing a merging event between two cold dark 
matter halos. Since the dark matter halos are much more 
extended than the bright cores of the galaxies, the interac- 
tion between the dark halos is already quite strong, while 
at the same time the stellar distribution is much less af- 
fected. However, by looking carefully at the outer part of 
the luminous halo's of the galaxies, some indications of in- 
teraction are visible. In particular, the lack of asymmetry 
in the outer contours of the light distribution indicates that 
the small galaxy at the center of the group has probably 
been stripped of its outer parts. This issue will be ana- 
lyzed in more detail in a forthcoming paper. The analysis 
of this particular lens makes it clear that the combination 
of the perturbative method and good HST data open inter- 
esting new possibilities. Using this approach it should be 
possible to probe the structure of the dark component at 
the mass scale of galaxies for a large number of systems. 
Statistical results about the geometry and structure of the 
dark halos should be derived, and will probably offer the 
possibility to evaluate the amount of substructure in the ha- 
los (Alard 2008). It is clear also that the mass-luminosity 
relation could be analyzed for a number of other systems, 
and that the discovery of systems similar to this one would 
constrain even more a modification of gravity. 



6.3. Density 

To estimate the density of the lens it is necessary to make 
some assumptions. Making a specific assumption about the 
background allows us to make an estimation of K2 and also 
break the mass-sheet degeneracy. For nearly isothermal lens 
distributions, the density iso-contours are given by Eq. ID. 21 
Note that in this equation the terms corresponding to the 
impact parameters cancel out 



dr,. 



d'fo ^ , d^fo 
+ /i — 



d^9 



d^9 



To optimize the signal to noise ratio, the iso-contours pre- 
sented in Fig. [10] are reconstructed by using an order 4 
Fourier expansion (see Sec. 4.6.3 for details). Note that a 
method to reconstruct the iso-contours of a nearly isother- 
mal lens was also studied by Evans & Witt (2003). However, 
their description of the potential is less general since it con- 
tains only one angular functional, it is suitable only for po- 
tentials with constant isophotes. 



Appendix A: Physical meaning of the perturbative 
approximation. 

For circular sources, Eq ( [8|) shows that /i is the mean 
radial position of the image. The width of the image 



7?2 



dfo 

de 



The condition for image formation is 



l^^l < i?o, thus images are formed near the minimas of 
|^|. The minimas of |^| are of two different types: lin- 
ear behavior of near the minima (small images), or a 
higher order polynomial behavior: caustics. To illustrate 
these properties, let's take the simple case of an elliptical 
potential, with ellipticity parameter rj: 

cj)^F + (1 + ^)2/') 



Let's now develop this potential to first order in 77: 
4> ~ F{r) - -F'{r)rcos29 



(A.l) 



7. Discussion 

This paper describes the structure of the dark matter enve- 
lope of a small group of galaxies. It is shown that the dis- 
tribution of the dark mass does not follow light (Fig. [T0|) . 
This is direct evidence for an independent dark component 
at the scale of galaxies. Such new results about a system 
with a mass much lower than a typical cluster (for instance 
the Bullet Cluster) are essential. First because there has 
always been a problem with missing baryonic mass in clus- 
ters (for instance X ray gas); second because structures like 
the Bullet Cluster contains many galaxies, and their dis- 
tances are not known with good enough accuracy to avoid 



As a consequence, (j)Q{r) ~ F{r), il}{r,9) = ^F'{r)rcos29 
(see Eq. \^ and since the background potential must be 
critical at r = 1, F'(l) = 1. Using Eqs ( [5|) and ( [7]) the 
perturbative fields reads: 



-Ml- 



2 V- , -F"(l))cos26l-t-xocos6'-|-?/osin6' 
7] sin 29 — Xq sin 9 + yo cos 9 

(A.2) 

Placing the source on the X-axis (yo = 0), at an impact 
parameter a;o — 2r/ Eq. f lA.2p one obtains at 9 = 0: — 

0, n = 1..3, which corresponds to a broad minima of |^^| 
near 9 = (see Fig. IA.1[) . This minimum corresponds to 




Fig. 10. The upper plot presents the potential iso-contours, the lower plot presents the density iso-contours. In the upper 
plot, the red curves represent the two 3 a contours of the potential iso-contours. The dashed curve is the unit circle. 
The blue curve represents the local deviation from circularity of the potential. The distance between the blue curve and 
the unit circle is proportional to the deviation from circularity. The lower plot presents the two 3 a iso-contours of the 
density. In each plot the small elliptical contours represent the galaxies in the lens. The small black contours represent 
the galaxies. 
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lens p by the following formula: 



T3 
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Fig. A.l. Plot of the function ^ [xq ~ 2rj and = in 

Eq. IA.2p . Images are formed when ^ < i?o- The nature of 

the image is related to the behavior of | \ near its minima, 
the arc corresponds to a large minima, with cancellation 

of the first and second order derivatives of For the 

small image, the behavior of near zero is linear, and 

corresponds to a much smaller minima of | | . 



■ihl Jo Io~ ^('"' sinnw u""*"^ du dv 



(B.2) 



Using Eqs ( IB.ip and ( IB.2p . and noting that: 

V /[r=l] V J[r=l] 

[oe) 



de 



[r=l] 



' /i = I]n " («« ^ c„) COS nO + 71 [hn - dn) siu n9 



J2n + '^n) '^'-"^ + ('^« + '^") '^^ 



(B.3) 

As a consequence, there is a direct relation between the 
Fourier expansion of the perturbative fields and the multi- 
pole expansion of the potential at the critical circle. This 
relation also demonstrates that the Fourier expansion of the 
perturbative fields is not equivalent to the Fourier expan- 
sion of the potential by Trotter et al. (2000). The expansion 
proposed by Trotter requires a number of additional coef- 
ficients with respect to the multipole expansion. On the 
other hand, Eqs dRT]) . {MM and (EH) indicate that the 
perturbative expansion is directly related to the multipole 
expansion without needing any additional coefficients. It is 
also clear from Eqs ( IB.ip and ( IB.3P that since the Fourier 
expansion of elliptical potentials is approximately of order 
2 (see Eq 



lA.ll) . the corresponding Fourier expansion will 
also be of order 2. Higher order terms in the Fourier expan- 
sion of the perturbative fields will thus be related to the 
skewness, boxyness, etc. , of the potential. The first order 
terms are related to the centering of the potential. 



the formation of an arc. This is a case of a cusp caustic, and 

behaves like an order 3 polynomial near the singularity. 
The counter-image in the opposite direction is much smaller 

and corresponds to a narrow minimum of |^^|- For this 

small image — for n — 1 only, and the behavior of 
this field is linear near the image. 

Appendix B: Relation to multipole expansion 

The expansion of the perturbative potential at a given 
radial position r ip reads (Kochanek 1991): 

, an(r) bn{r) . 
w = — > COS nO H sm nO+ 

Cn(?') r"^ COS nO + dn{r) r"sinn6' 

The perturbative theory requires the evaluation of the po- 
tential gradient on the critical circle. Since the coordinate 
system has been re-normalized so that the critical radius is 
situated at r = 1, all the evaluation of quantities related to 
the potential will be performed at r = 1. The coefficients 
(a„, 6„, c„, (i„) are related to the projected density of the 



Appendix C: Geometrical interpretation 

The perturbative fields are directly related to the geomet- 
rical structure of the potential. The potential iso-contour 
near the critical circle is easy to infer from Eq. ( ID.ip . The 
iso-contour equation reads: 

Mr)+eh{e)+eh{e){r~l)^C 

Searching for the iso-contour near the critical circle = 
1 -|- edri and developing the former equation to first order 
in e: 

e{dn + ,h)=C -Ml) 

When no perturbation is introduced (/o = 0), the iso- 
contour must be reduced to the critical circle {dr = 0), 
as a consequence the right side of the equation must be 
zero too, and the former iso-contour equation reduces to: 

dn - -/o (C.l) 

The potential iso-contour depend only on /o, while the iso- 
contour variation depends also on /i. Considering a con- 
tour with radial position r^, defined by (f>{ri,9) = C, the 
equation of a nearby contour is cj)(ri + Sri, 9) — C + dC, 
or equivalently ^{ri,9)Sri — dC. There is no iso-contour 
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variation provided that Sri is constant, which is equivalent 
to ^{ri,9) being constant. To first order in e: 



deb 



1 + /i - 



-(l)/o 



(C.2) 



Considering that Eq ( IC.2|) depends only on /i for purely 
isothermal backgrounds, and that CDM halos are quite 
close to isothermal in general, the iso-contour variation is 
dominated by /i. 

Appendix D: Relation to the geometry of the 
density 

The two perturbative fields are related to the potential by 
the following expansion (Eq. d]): 



</> = (/.o(r)+e/o(0) + e/i(0)(r-l) 



(D.l) 



Note that since r — 1 + edr the expansion presented in Eq. 
( ID.1|1 is of order 2 in e. This expansion reduces to order 
1 when introduced in the lens equation. By inserting Eq 
( ID.ip in the Poisson equation, taking r = 1 + edr, and 
developing to first order in e, the density reads 



p ~ (1 + 2C2) + (2^2 + 6C3 - l)dr + 



d^e 



+ h +2/2 6 



The iso-contour near the critical circle is defined by the 
equation p = C. In the absence of a perturbation {fi = 0) 
the iso-contour must reduce to the critical circle {dr = 0). 
As a consequence the equation of the iso-contour is: 



(1 - 2C2 - 6C3)dr 



d^d 



+ h + 2/2 



For a nearly isothermal distribution, the higher order 
derivatives (n > 1) are close to zero. Assuming that these 
higher order derivatives are of order e (C2 — C3 ~ e and 
/2 ~ eF{9)), the former equation reduce to: 



dri 



(D.2) 
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